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I. INTRODUCTION 

Besides classical and semi-classical descriptions of dissipative molecular systems several 
quantum theories exist which fully account for the quantum effects in dissipative dynamics. 
Among the latter are the reduced density matrix (RDM) formalism |jl]-§l and the path 
integral methods |P,0]. Here we concentrate on the Redfield theory 0,0, in which one 
has to solve a master equation for the RDM. It is obtained by performing a second order 
perturbation treatment in the system-bath coupling as well as restricting the calculation to 
the Markovian limit. With this approach the quantum dynamics of an "open" system, e.g., 
the exchange of energy and phases with the surroundings modeled as a heat bath, can be 
described. The unidirectional energy flow into the environment is called dissipation. Within 
this theory it is possible e.g. to simulate the dissipative short-time population dynamics 
usually detected by modern ultrafast spectroscopies. 

Because the Redfield theory is a Markovian theory the time evolution of the RDM is 
governed by equations containing no memory kernel. In the original Redfield theory |p,|^ 
the secular approximation was performed. In this approximation it is assumed that every 
element of the RDM in the energy representation is coupled only to those elements that 
oscillate at the same frequency. In the present study we do not perform this additional 
approximation. 

For larger systems numerical implementations for solving the Redfield equation are nu- 
merically very demanding and therefore require that one finds the most appropriate way 
to perform the time evolution of the RDM. Straightforward one can construct and directly 
diagonalize the Liouville superoperator. For processes with a time-independent Hamilto- 
nian the rates, i.e. the characteristic inverse times of an exponential decay of the occupation 
probability of the excited states, can be obtained in this way. In such an approach a huge 
number of floating point operations will be involved and the overall computational effort 
will scale as M^ where M is the size of the basis of vibronic states. Furthermore the direct 
diagonalization can be numerically unstable, but nevertheless has been successfully used (see 
e.g. fl^])- Another strategy suggests solving N"^ ordinary differential equations and requires 



products between the Liouville superoperator and the RDM which scale as Af^ (see e.g. [0) 
This is also numerically demanding for larger systems. Assuming a bilinear system-bath cou- 
pling the numerical effort can be reduced considerably by rewriting the Redfield equation 



in such a form that only matrix-matrix multiplications are needed |T2| rather than applying 
a superoperator onto the RDM. Hence, a computational time scaling of M^ and a storage 
requirement of N"^ is achieved. In the present paper our numerical studies will be based on 
this approach. 

If one wants to reduce the scaling of the numerical effort with increasing number of 
basis functions even more one has to go to stochastic wave function methods ||13|-Q . They 



prescribe certain recipes to unravel the Redfield equation and to substitute the RDM by a 
set of wave functions which evolve partially stochastically in time. The method will have 
the typical scaling of the well developed and optimized wave function propagators, i.e. N"^. 
It has been applied, for example, to electron transfer systems [p!^ -p!9[| and shown to give 



accurate results. Direct solutions and the stochastic wave packet simulations have already 
been compared numerically ||2y, ^]. All these first studies were restricted to dissipation 



operators with Lindblad form [^. Breuer et al. |2^ showed that the stochastic wave function 



approach can also be applied to Redfield master equations without the secular approximation 
and for non-Markovian quantum dissipation. In particular for complex systems with a 
large number of levels its practical application is very advantageous. Between the direct 
and the stochastic methods to solve the Redfield equation the accuracy differs especially 
because direct RDM integrators are numerically "exact" while the stochastic wave function 
simulation methods have statistical error. The scope of this work are small and medium size 
problems. Therefore we compare the different numerical algorithms for a direct integration 
of the Redfield equation. But when using stochastic methods for density matrix propagation 
one has to solve Schrodinger-type equations with a non-Hermitian Hamiltonian. For that 
purpose the same algorithms as investigated here can be used. In this sense the present 
study is also of importance for solving Redfield equations by means of stochastic methods. 
Here we use different numerical schemes to solve the Liouville-von Neumann equation. 
The performance of the well-known Runge-Kutta (RK) scheme is studied in two different 
implementations: as given in the Numerical Recipes [^ and by the Numerical Algorithms 
Group 1^. Compared to these general-purpose solvers are more special algorithms which 



have been applied previously to the time evolution of wave packets and density matrices. For 
density matrices these are the Short-Iterative- Arnoldi (SIA) propagator [|1^,^, the short- 
time Chebyshev polynomial (CP) propagator |TT|, and the Newtonian polynomial (NP) 



propagator ^,|2^. The latter propagator is also used as a reference method because of 
its high accuracy. In addition, a symplectic integrator (SI), which was originally developed 
for solving classical equations of motion and extended to wave packet p9|,p0[ and density 



matrix propagation |]31[, is tested. 

Besides the propagation algorithm one has to determine the appropriate representation 
in which to calculate the elements of the RDM and the Liouville superoperator. The choice 
depends strongly on the type of physical problem that one considers. Coordinate (grid) rep- 
resentation has an advantage when dealing with complicated potentials, e.g. non-bonding 
potentials. For the problem of dissociation dynamics in the condensed phase the grid repre- 
sentation has been applied based on a Lindblad-type master equation [^]. Other examples 



using grids are works of Gao 1^31 and Berman et al. |^ . Convenient treatment of electron 



transfer dynamics is done in state representation, because one can model the system using 



a set of harmonic diabatic potentials [^. Other authors [jT2|,^^ choose the adiabatic 
eigenstates of the whole system as a basis set to treat similar problems. A comparative anal- 
ysis of the benefits and drawbacks of the diabatic and adiabatic representation in Redfield 



theory will be given elsewhere [^|. From the viewpoint of numerical efficiency we focus on 
both representations in the present article. 

The paper is organized as follows. In the next section we will make a short introduction 
to the model system and a discussion on the versions of the Redfield equation and the 



numerical scaling that they exhibit. In Section |IT| the methods for propagation used in this 
work are briefly reviewed. In Section ^ we compare the efficiency of several propagators in 
solving a simple electron transfer problem with multiple levels. A summary is given in the 
last section. 



II. THE REDFIELD EQUATION AND ITS SCALING PROPERTIES 

In the RDM theory the full system is divided into a relevant system part and a heat 
bath. Therefore the total Hamiltonian consists of three terms - the system part Hs, the 
bath part H^ and the system-bath interaction -f^sB^ 

H = Hs + Hb + Hsb. (2.1) 

The separation into system and bath allows one to formulate the system dynamics, 
described by the Liouville-von-Neumann equation, in terms of the degrees of freedom of 
the relevant system. In this way one loses the mostly unimportant knowledge of the bath 
dynamics but gains a great reduction in the size of the problem. Such a reduction together 
with a second order perturbative treatment of the system-bath interaction H^b and the 
Markov approximation leads to the Redfield equation |jl],H,|,||: 

p=-^[Hs,p]+np = Cp. (2.2) 

In this equation p denotes the reduced density operator and TZ the Redfield tensor. If one 
assumes bilinear system-bath coupling with system part K and bath part $ 

^SB = K^ (2.3) 

one can take advantage of the following decomposition [Ellir 



p=~ [Hs,p] + ^{[Ap,K] + [K,pA^]}. (2.4) 

Here the system part K of the system-bath interaction and A together hold equivalent 
information as the Redfield tensor TZ. The A operator can be written in the form 

A = K f dT{^{T)^{0))e'^^ = KC{lu) . (2.5) 



The two-time correlation function of the bath operator C(r) = ($(r)$(0)) and its 
Fourier-Laplace image C{uj) can be relatively arbitrarily defined and depend on a micro- 
scopic model of the environment. Different classical and quantum bath models exist. Here 



we take a quantum bath |2^, i.e. a large collection of harmonic oscillators in equilibrium. 



that is characterized by a Bose-Einstein distribution and a spectral density J{uj): 



C{uj) = 27r 



1 + (e^'^/^^^ - l) ' 



[j(^) _ j(_^)] . (2.6) 



Now we look for bases that span the degrees of freedom of the relevant system. Let us 
consider atomic or molecular centers m at which the electronic states \m) of the system 
are localized. Their potential energy surfaces (PESs) will be approximated by harmonic 
oscillator potentials, displaced along the reaction coordinate q of the system. As such a 
coordinate one can, for example, choose a normal mode of the relevant system part which 
is supposed to be strongly coupled to the electronic states. The centers are coupled to each 



other with constant couphng Vmn- For example two such coupled centers are sketched in 
Fig. |l|. The coupled surfaces |1) and |2) are assumed to describe excited electronic states. 
The electron transfer takes place after an excitation of the system from its ground state |g). 
Using this microscopic concept we define K as the system's coordinate operator 



K = q = J2i'^^mM/n) 



m ' ^iTi 



|r7i)(r?7,| 



(2.7) 



where am and a\^ are the boson operators for the normal modes at the center m, M. is the 
reduced mass and ujm are the eigenfrequencies of the oscillators. In the same picture the 
Hamiltonian of the relevant system reads 



H^ = Ei'^" 



Um + ( aL«m + 2 ) ^^rn 



1 



Vmn]\m){n\ . 



(2.8) 



From this point on we consider two possible state representations in order to calculate the 
matrix elements of p and the operators ifs? K and A. The diabatic (local) basis is a direct 
product of the eigenstates \M) of the harmonic oscillators and the relevant electronic states 
\m) (Fig. m, left panel). The intercenter coupling Vmn gives rise to off-diagonal elements of 
the Hamiltonian matrix 



{mM\H^\nN) 



U„ 



M +-\huOr. 



^mn^MN + (1 " ^mn) Vmn{'mM\nN) 



(2.9) 



where Um are the energies of the minima of the diabatic PESs. Other important properties 
of the diabatic representation are the equidistance in the level structure and the diagonal 
form of the system-bath interaction operator H^b- This determines the tridiagonal band 
form of the operator K. 

When neglecting the influence of the intercenter coupling on the dissipation a very effi- 
cient numerical algorithm can be derived [p4| , |38| . Of course, for strong intercenter couplings 
the populations {mM\p\mM) at long times deviate from their expected equilibrium values. 
But even for very small couplings there are cases in which the population does not converge 
to its equilibrium value [^ . Therefore this neglect of the influence of the intercenter cou- 
pling on the dissipation has to be handled with care. On the other hand this approximation 
makes the extension of the present electron transfer model to many modes conceptually 



much easier 11171,^. 



The matrix elements of the operators in the dissipative part of Eq. ( p.4|) read 
{mM\K\nN) = {2uJmM/ny^'^ 5mn (6m+i,nVM + 1 + 6m-i,nVm) 
and 

{mM\A\nN) = {mM\K\nN)C{umMnN) ■ 
In Eq. ( |2.11| ) uj^muN denote the transition frequencies of the system 

humMnN = {mM\Hs\mM) - {nN\Hs\nN). 



(2.10) 



(2.11) 



(2.12) 



Since the system can emit or absorb only at the eigenfrequencies of the system oscillators 
Um the spectral density of the bath J{uj) is effectively reduced to a few discrete values 



J{uj) = J2mlmS{^ — ^m)- The advantage of this approach hes in the scahng behavior with 
the number of basis functions. As shown in Fig. ^ it scales hke Af^'^ where A/" is the number 
of basis functions. This is far better than the scahng without neglecting the influence of the 
intercenter coupling on the dissipation as described below. 

If the system Hamiltonian Hs is diagonalized it is possible to use its eigenstates as a basis 
(Fig. ID, right panel), in which to calculate the elements of the operators in Eq. (|2.4| ). Of 
course there will be no longer any convenient structure in K or A, so that the full matrix- 
matrix multiplications are inevitable. For this reason the computation of Cp{t) scales as 
A/"^, where A/" is the number of eigenstates of H^. There appears to be a minimal number 
Ao below which the diagonalization of H^ fails or the completeness relation for \mM) is 
violated. Nevertheless, the benefit of this choice is the exact treatment of the system-bath 
interaction. Denoting the unitary transformation that diagonalizes Hs by U one has 

e = U^HsU (2.13) 

where e is a diagonal matrix containing the eigenvalues. In this way it is straightforward to 
obtain the matrices for p and K. Equation ( p.6|) for C{uj) still holds but a new definition of 
the spectral density J{uj) is necessary because of the non-equidistant adiabatic eigenstates. 
The bath absorbs over a large region of frequencies and this is characterized in the model 
by J{uj). One needs the full frequency dependence of J{uj) which we take to be of Ohmic 
form with exponential cut-off: 

J{uj) = r]e{uj)uje-'^/''\ (2.14) 

Here 6 denotes the step function and Uc the cut-off frequency. In this study all system 
oscillators have the same frequency Ui (see Table |) and the cut-off frequency Uc is set to be 
equal to cui. The normalization prefactor t] is determined such that 

diuJ{uj) = -fi . (2.15) 



Equation ( p.l4|) together with Eq. (|2.6| ) yields the correlation function in adiabatic repre- 
sentation. 

The introduced representations allow us to consider the numerical effort for a single 
computation of the right hand side of Eq. (|2.4|), i.e. of Cp(t). In the diabatic representation 



its computation can be approached using two different algorithms. It is possible to per- 
form matrix-matrix multiplication only on those elements of K and A which have nonzero 
contributions to the elements of Cp{t) (Fig. ^ solid line). This is advantageous because of 
the tridiagonal form of K in diabatic representation and shows the best scaling properties, 
namely Af^'^. In the same representation but performing the full matrix-matrix multiplica- 
tions in Eq. ( p.4|) (Fig. H, dashed line) the scaling behavior is slightly worse than the same 
operation in adiabatic representation (Fig. |^, dotted line). This is due to the non-diagonal 
Hamiltonian Hs in the former case that makes the computing of the coherent term (see 
Eq. ( p.2| )) more expensive. Below we will take the full matrix-matrix multiplication to eval- 
uate Cp{t) in both representations. We do this to concentrate on the various propagation 
schemes, not the unequal representations. Nevertheless there are performance changes in the 
different representations because of the disparate basis functions and forms of the operators 
in these basis functions. 
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III. THE DIFFERENT PROPAGATION SCHEMES 

A. Runge-Kutta method 

The RK algorithm is a well-known tool for solving ordinary differential equations. Thus, 
this method can be successfully applied to solve a set of ordinary differential equations for 
the matrix elements of Eq. ( p.4|) . It is based on a few terms of the Taylor series expansion. 
In the present work we use the FORTRAN?? implementation as given in the Numerical 
Recipes p^ which is a fifth-order Cash-Karp RK algorithm and will be denoted as RK-NR. 



As alternative the RK subroutine from the Numerical Algorithms Group [EH] which is based 



on RKSUITE |3|] was tested with a 4(5) pair. It will be referred to as RK-NAG. Both RK- 
NAG and RK-NR involve terms of fifth order and use a prespecified tolerance r as an input 
parameter for the time step control. The tolerance r and the accuracy of the calculation are 
not always simply proportional. Usually decreasing r results in longer CPU times. 

In our previous work |i^ a time step control mechanism different from those used in RK- 



NAG and RK-NR was tested. Discretizing the time derivative in Eq. ( p.2|) and requiring 

p{U+i) - p{ti 



At 



+ Cp{ti) 



< T (3.1) 



one only has to call the propagation subroutine once and to store the previous RDM. In 
addition one has to calculate the action of the Liouville superoperator L onto the RDM but 
the numerical effort for this is small compared to a call of the propagation subroutine. It was 
shown in Ref. |^0[ that this time step control is the most efficient for propagation with the 
coherent terms in Eq. (|2.2| ) only but disadvantageous for problems with dissipation. This is 
the reason why we do not include this algorithm in the present study. 

B. Short Iterative Arnold! propagator 



The SIA propagator ||T2|,^ is a generalized version of the Short Iterative Lanczos propa- 
gator |^I| to non-Hermitian operators. With the Short Iterative Lanczos algorithm the wave 
function can be propagated by approximating the time evolution operator in Krylov space, 
which is generated by consecutive multiplications of the Hamiltonian on the wave function. 
In analogy the Krylov space within the SIA method is constructed by recursive applications 
of the Liouville superoperator onto the RDM p„ = £"p(t). In this way it is tailored for the 
RDM at every moment in time. The Liouville superoperator, denoted by I in Krylov space, 
has Hessenberg form 

L ^ VIV^ , (3.2) 

where the orthogonal transformation matrix V is constructed iteratively using the so-called 
Lanczos procedure WM- The Krylov representation / can be easily diagonalized to L with 



the help of a transformation matrix S: 

e^* ^ VSe'^'S-W^ . (3.3) 

Since the diagonalization is performed in the Krylov space the numerical effort depends 
on its dimension which can be chosen small in practice. Having thus derived a diagonal 
operator e^^ the calculation of p{t) is straightforward. 



C. Symplectic integrator 



The Sis were originally developed for solving classical equations of motion |^2[. The time 
evolution of a classical Hamiltonian system can be viewed as a canonical transformation 
and Sis are sequences of canonical transformations. Recently it was shown that the time 
evolution of wave packets [^,|3y] and density matrices |3T| can also be performed using Sis. 
In order to rewrite the Redfield equations in the form of coupled canonical variables that 
are analogous to classical equations of motion one defines the functions [BT 



Q{t) = p{t) 

Pit) = m 



the operator 



and the Hamiltonian function 



W 



1 



-—£2 

-.2'^ 



(3.4) 
(3.5) 



(3.6) 



G{Q,P) = -[P^P + Q^WQ]. 
Doing so one obtains equations of motion analogous to the classical ones 



d 
di 
d 



Pit) 



dGiQ,P) 
dQ 



-WQ{t) 



Rewriting this into the SI algorithm of order m yields 

6,- At 



Pr = P^^l + 



Qi — Qi-1 + 



3 *- ^i — 1 



h 
ttiAt 

IT 



(3.7) 

(3.8) 
(3.9) 

(3.10) 
(3.11) 



for i = 1, . . . ,m. Different sets of coefficients {oj} and {bi} are given in the literature. Here 
we choose the McLachlan-Atela fourth-order method [E3l. The coefficients for this method 



are listed in Ref. |^^. A comparison of the McLachlan-Atela fourth-order method with the 
McLachlan-Atela third-order method p3| and Ruth's third-order method P2[ has been given 



elsewhere 31 . 



D. Newton polynomial scheme 



Another way to solve Eq. ( |2.2| ) is by a polynomial expansion of the time-evolution opera- 
tor. Such methods are well established and approved for wave-function propagation |^,41 . 
Recently the Faber [^ and NP [^ algorithms have been applied to propagate density ma- 
trices and it has been shown that they behave very similarly EH] . The main idea of the NP 



method is the representation of the Liouville superoperator by a polynomial interpolation 



n-1 



M 



e- ^ Vn,-i{^) = E «nP„ = E «n IK-^ - ^j) (3-12) 

n=0 n=0 j=0 

of order Np where the p„ are computed recursively and a„ are the n-th divided differences. 
The interpolation points Xj can be chosen to form a rectangular area in the complex plane 
(see Fig. ^ which contains all eigenvalues oi C. This interpolation scheme is uniform, i.e., the 
accuracy in energy space is approximately the same in the whole spectral range of C. This is 
in contrast to schemes such as the SIA propagator which are nonuniform approximations. A 
consequence of this property is the very high accuracy which can be achieved with uniform 
propagators. This is why we take a high-order NP expansion as reference solution. Since the 
quality of the approximation of the time evolution operator is equivalent to a scalar function 
with the same interpolation points Xj, one can, before performing the actual calculation, 
check the accuracy on a scalar function. For the calculation with the NP propagator we set 
the truncation limit of the expansion to 10~^^, i.e., the sum in Eq. ( |3.12| ) is truncated when 



the residuum fulfills ctrillPnll < 10 [^ 



E. Chebyshev polynomial scheme 

As a last contribution to the present study we will examine the CP propagator. Re- 
cently it was studied by Guo et al. [|ll[] for density matrices. The Liouville superoperator is 



approximated by a series of CPs Tk{x). Generally the CPs diverge for non-real arguments. 
For propagators of the kind e"*''^* it has been shown ^8\ that the CPs may tolerate some 



imaginary part in the eigenvalues of H. The stability region has the form of an ellipse 
with a center at the origin and a very small half-axis in imaginary directions p^. In con- 



trast, the eigenvalues of the Liouville superoperator are spread over the negative real half 
of the complex plane and symmetrically with respect to the real axis (see Fig. |^). The real 
components for the system that we consider are one order of magnitude smaller than the 
imaginary components. This is why we make the expansion along the imaginary axis and 
use an expression similar to that already applied to wave function propagation pT|: 



^ct ^ ^L+At ^ ^2 - Sr.o)UL~At)Tn{L) . (3.13) 

n=0 

Here the expansion coefficients J„ are the Bessel functions of the first kind, and L is the 
appropriately scaled Liouville superoperator: L = {C — L^)/L~ , where L~ and L~^ are the 
half span and the middle point of the spectrum of C. Since the spectrum is symmetric with 
respect to the real axis, L~^ = 0. The time evolution of p is given by 

Np-l 

pit + At) ^ E (2 - 5no)Jn{L-At)prr . (3.14) 

n=0 

The Chebyshev vectors pn are generated by means of a recurrence procedure: 

Pn = 2Lpn-i + Pn-2, Po = p{t) and Pi = LpQ . (3.15) 
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For the CP and NP methods one has to adjust the values of the spectral parameters L~ 
and L"*". One can obtain some knowledge about the spectrum of C by an approximate diag- 
onalization, e.g. by Krylov subspace methods. For instance, Fig. ^ shows an approximate 
spectrum of £, appropriately scaled so that all eigenvalues lie within the rectangle formed 
by the Newtonian interpolation points. 

IV. PERFORMANCE OF PROPAGATION METHODS 

The aim of this section is to compare the different numerical methods described above 
for propagating the RDM in time. The calculations were performed for both RK methods 
with different tolerance parameters r and for the SI as well as the NP, CP, SIA propagators 
with different timesteps. The number of expansion terms Np in NP and CP propagators is 
170 and 64, respectively. The dimension of the Krylov space for the SIA method was set 
to 12 because smaller as well as larger values are less efficient for the example studied here. 
All computations were made on Pentium III 550 MHz personal computers with intensive 
use of BLAS and LAPACK libraries. The code was compiled using the PGF90 Fortran 
compiler |^ . For estimation of the computational error of all methods the NP algorithm 



with 210 terms was chosen as a benchmark. 

In this work we consider only two centers m = 1, 2 which is the minimal model to describe 
the main physics of an electron transfer reaction. A basis size of 16 levels per center satisfies 
the completeness relation and presents no difficulties during the diagonalization of H^- The 
electronic couphng was fi2 = 0.1 eV. We choose 7i = 72 = 1.57863 x 10~^ eVA~^ and 
M. = 20mp where rrip is the proton mass. The temperature T = 298 K is used. In Table | 
the parameters for the system oscillators are given. 

The process that is simulated involves the following scenario. A Gaussian wave packet 
is prepared as initial state by a vertical transition from the lowest vibrational level of the 
ground electronic state |g) to the first (upper) center |1): 

PiMiNit = 0) = (lM|gO)(gO|liV) . (4.1) 

The energy distribution of the occupied eigenstates by the wave packet depends on the 
displacement q^ — q^ between the PESs of |g) and |1). During the pulse the two excited 
electronic states |1) and |2) are assumed to be decoupled. In this way one can simulate 
the absorption of electromagnetic radiation from a pulse with vanishing width. Right after 
the pulse is over, the wave packet starts moving on the excited PESs and spreading. The 
relevant system part begins losing energy to the bath and dephasing. The population on 
the upper center starts decaying. When the damping is not too strong, as for the model 
parameter studied here, a damped oscillation of the population between the two excited 
PESs can be seen. We assume no coupling to the ground state after the pulse. After a 
certain time the system reaches its equilibrium state. 

In all cases the RDM was propagated for a total time period of 3 x 10^ a.u. which is 
sufficient for complete relaxation to equilibrium. It was compared to the RDM p^ef evaluated 
by the NP algorithm at the same points in time. The relative error e{t) of each method at 
a certain moment in time t has been estimated using a formula similar to that proposed for 
wave functions by Leforestier et al. | |41[ |: 
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e{t) 



^_TT{p{t)p,,,{t)) 



Tr(pLW) 



(4.2) 



As the error e we define the maximum value oie{t) over the total propagation time. For more 
details we refer to our previous paper |40|. Other error measures (see for example [11, 47|) 



can be used as well but they will have the same qualitative behavior. 

As an index for the numerical effort two possibilities were explored. The first one is a 
direct measurement of the CPU time of the total propagation (Fig. ^. It may look quite 
different on other computer architectures or even on the same architecture but under changed 
operation conditions. An evidence of the performance (Fig. 1) will be expressed by means 
of CPU time versus the error e. 

Another approach to describe the numerical effort has been proposed [|l^ and is called 



efficiency factor. It is defined as the ratio between the timestep At and the number of 
operations Cp{t) within this timestep. Because of the definition it is a machine independent 
quantity. The larger the efficiency factor, the better the performance of the algorithm. 
Because the RK algorithms propagate with variable timestep we cannot directly use the 
definition of the efficiency factor. Instead we define a quantity a as the total number Nc of 
£p(t)-evaluations divided by the total time for the propagation: 

a = Ar,/(iV,At) . (4.3) 

Here Ng denotes the total number of timesteps. The inverse of a will have the meaning of 
an efficiency factor for an averaged timestep At. We should point out that N^ does not take 
into account the effort for summation of the different contributions. In particular in the 
case of the NP method the summation of the different terms in the polynomial expansion 
Eq. ( |3.12| ) can be non-negligible. This can be seen in the different relative performance of 
the propagators shown in Figs. ^ and p. We consider both the CPU time and the quantity 
a as measures of the numerical effort. 

Contributions from the algorithm to calculate Cp{t) also influence the CPU time. As 
discussed above, in all computations represented in Figs. § and ^ the full matrix-matrix 



multiplications in Eq. (|2.4|) were performed. The performance of the CP, NP, SI and SIA 
methods is only influenced very little by the choice between diabatic and adiabatic repre- 
sentation. Both RK implementations are less efficient in the adiabatic than in the diabatic 
representation, though the RK-NAG scheme has still the best performance besides the NP 
algorithm. The RK-NR scheme has an advantage for computation in diabatic rather than 
in adiabatic representation especially for medium precision requirements. In that range the 
performance curves of the RK methods exhibit a shoulder for the adiabatic case which seems 
to result from a numerical artifact. 



Because the error of the SIA algorithm is not uniformly distributed in energy space |^ 
we could expect some difference in its performance in diabatic and adiabatic representation. 
But because the coupling v^ chosen here is not very large, the eigenstates of the coupled 
system are just slightly disordered (see Fig. |^, right plot) and hence the performance of the 
SIA algorithm is almost not changed. 

The uniformity, stability and high accuracy of the CP propagator for wave functions is 
well known |^, ^, Q . The CP approach to density matrix propagation was introduced 
by Guo and Chen [^. Using a damped harmonic oscillator as model system and starting 
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from a pure state they established that the relative error can reach the machine precision 
limits (10~^^) for sufficiently small stepsize. However, for the system of coupled harmonic 
potentials studied here and using an initial RDM with non-zero off-diagonal elements the 
error saturates at e ~ 10~^ (see Fig. ^). It was not possible to decrease this saturation 
limit of e neither by increasing the order of the CP nor by decreasing the timestep. This 
saturation limit seems to depend strongly on the imaginary part of the eigenvalues of C For 
large timesteps the CP method loses its stability and one needs to estimate the efficiency 
range of Np, At and L~. Turning off the dissipation we could reach much higher accuracy 
with the CP propagator as expected. 

The SI is easy to implement. The expansion coefficients are fixed and can be taken 
from literature. At the same time the fixed coefficients seem to limit the accuracy. For not 
too high accuracy the performance of the SI is as good as that of the other propagators in 
adiabatic representation. In diabatic representation its performance is a little worse. But 
we were not able to achieve very high accuracy with this method. This might be due to the 
special version, the fourth order McLachlan-Atela method, which we chose. 



As already highlighted ^^ the NP scheme is very stable for arbitrary spectral properties 
of jC. The only restriction is that the spectrum must be confined within the area formed 
by the interpolation points. In our investigation the NP propagator performs with a good 
accuracy for timesteps of 1500 a.u. {Np = 170) which is 10 times larger than the step size 
of the CP scheme. Higher order expansions might be even more efficient but the numerical 
implementation gets tricky and easily unstable. For timesteps of 100 a.u. and A''^ = 50 
the NP algorithm is already numerically exact but computational very expensive (see the 
arrows in Fig. ^ and Fig. ^. For problems with time-dependent Hamiltonians (e.g. non- 
stationary external fields with relatively small amplitude) the RK and SIA methods will be 
more efficient with small timestep. 

At the end we should point out that there exists no ultimate method to determine the 
performance of a certain numerical approach which could be valid for different platforms. 
Tuning and optimization features are generally not portable and this may cause even different 
scaling behavior and hence a different method of preference. That is the reason why the 
generality of the results is limited to similar computation platforms and even to systems 
with similar properties of the corresponding Liouville superoperator. But on the other hand 
this study can give hints on the performance of the different algorithms in general. 

V. SUMMARY 

In the present work an estimation of the numerical efficiency of several methods for 
density matrix propagation has been given. The example of electron transfer in a two-center 
system has been used for this purpose. A specific measure of the numerical effort has been 
introduced in order to compare methods with fixed timestep and such ones with timestep 
control (RK). Besides the method of reference (NP) the RK-NAG approach shows best 
performance for both cases of adiabatic and diabatic representation. The advantage of the 
SIA propagator is that the accuracy improves with decreasing the timestep in all cases we 
investigated. That is not the case with the CP propagator which exhibits a saturation of 
accuracy and is therefore not convenient for very small timesteps. The easy-to-implement SI 
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gives reasonable performance for not too high accuracy. The present SI seems to be hmited 
in accuracy due to the fixed coefficients. 

The present studies were restricted to state bases. Of course similar calculations can be 
done on a grid which is especially useful for complicated or unbound potentials. In these 
cases another propagator, the split operator |^, should be taken into account. This operator 
has the advantage that its performance does not (directly) depend on the spectral range of 
the Hamiltonian or Liouville operator. So it may perform very well for problems with a 
large spectral range although it cannot be applied to operators which have mixed terms 
in coordinate and momentum operators. The use of the mapped Fourier method [Q may 
reduce the number of grid points significantly and first wave packet propagations with this 



method have been done |51|,^. Recently the mult i- configuration time-dependent Hartree 
method has been established to treat density matrix operators [Q. This method might be 
favorable for multi-dimensional systems. 

The presented methods can be used for various applications in the field of dissipative 

molecular dynamics in condensed phases, where the RDM approach provides a good way 

of describing processes in systems with one or more degrees of freedom. This includes the 

electron transfer processes mentioned in the introduction as well as exciton transfer processes 

5^. It can also be used to simulate pump-probe experiments |TP|, surface scattering of atoms 



5^1 , etc. Also coherent control schemes in dissipative environments can be studied ||5^. So 
the numerical studies given here can be applied to a broad range of problems in physics, 
chemistry, and biology. 
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FIGURES 





FIG. 1. Potential energy surfaces (PESs) for a model electron transfer system. Diabatic PESs 
are plotted on the left side, and the PES of the adiabatic excited state |e) on the right side. 
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FIG. 2. Scaling behavior of the product Cp{t). Solid line - tridiagonal form of K in diabatic 
representation, dotted line - adiabatic representation, dashed line - diabatic representation with 
full matrix-matrix multiplications. The CPU time is scaled so that it is equal to 1 for AA = 50 in 
diabatic representation. 
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FIG. 3. Scaled spectrum L/L^ of the Liouville superoperator for the model of electron transfer. 
Approximate eigenvalues obtained in Krylov subspace are plotted as dots. Open squares denote 
the interpolation points Xj for the NP scheme. 
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FIG. 4. Numerical performance of different numerical propagators. Results obtained in the 
diabatic (adiabatic) representation are shown on the left (right) plot. The arrows represent the 
numerical performance for the NP propagator with 50 terms and timestep 100 a.u. 
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FIG. 5. Numerical performance of different numerical propagators. The numerical effort a is 



defined in Section IV . Results obtained in the diabatic (adiabatic) representation are shown on the 
left (right) panel. The arrows represent the numerical performance for the NP propagator with 50 
terms and timestep 100 a.u. 
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